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Abstract 

We introduce a scheme to describe the evolution of an interacting system of 
bosons, for which the field operator expansion is truncated after a finite number 
of modes, in a rigourously controlled manner. Using McLachlan's principle of 
least error, we find a self-consistent set of equations for the many-body state. 
As a particular benefit, and in distinction to previously proposed approaches, 
our approach allows for the dynamical increase of the number of orbitals during 
the temporal evolution. The additional orbitals, determined by the condition 
of least error of the truncated evolution relative to the exact one, are obtained 
from an initial trial state by a method we call steepest constrained descent. 
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1. Introduction 

When we treat an atom, such as 87 Rb, 23 Na, and 7 Li which consists of an 
even number of fermions, as a whole, the Fock space creation and annihilation 
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operators of these composites satisfy bosonic commutation relations. Since the 
experimental realization of Bose-Einstein condensates [l[ , a large variety of ex- 
periments with these bosonic isotopes have opened up a fascinating mesoscopic 
and macroscopic quantum world [2[. After an initial period, concentrating on 
the effective single-particle physics of these ultracold dilute gases Q , more re- 
cently their many-body physics came into focus, revealing rich and hitherto 
unexpected possibilities to test fundamental correlation properties at the mi- 
croscopic level Q. 

The extension to a true many-body physics, incorporating quantum cor- 
relations beyond mean-field, requires, however, vast computational resources 
when both the number of particles and the interaction between those particles 
increases. Therefore, a simplification of the problem by truncating the field 
operator expansion to a finite number of modes (or, as an equivalent term, 
single-particle orbitals) has been commonly utilized to obtain results relevant 
to the prediction of experiments in trapped bosonic quantum gases. The most 
extreme truncation, the semiclassical form of mean- field theory, retaining just 
one orbital, gives the well-known Gross-Pitaevskii equation. Without the aid of 
contemporary computers, it seemed to be inevitable until most recently, partic- 
ularly in out-of-equilibrium situations far away from the ground state, to reduce 
the complexity of the problem at hand as much as possible, and hence to use the 
Gross-Pitaevskii equation approach. With the increased interest in many-body 
physics, however, there arose the necessity to go beyond the all-too-simplificd 
mean field approach of the Gross-Pitaevskii equation. The accuracy of predic- 
tions on many-body correlations and the corresponding response functions will 
obviously increase with a less severe degree of truncation, though the solutions 
will not be exact still. 

To derive the equations of many-body evolution, various variational ap- 
proaches can be employed. Historically the first was the variational ansatz of 
Dirac and Frenkcl [111, followed by McLachlan's variational principle [7[ and the 
time-dependent variational principle (TDVP), which is a principle of stationary 
action [1,0. Therefore, there are various, not necessarily equivalent, choices of 
variational principle for finding the equation of motion of the truncated many- 
body evolution. The Dirac-Frenkel principle imposes (5&\H— idt\&) = (h = 1), 
where (5$ denotes any possible variations of the many-body state $ with respect 
to a given set of variational parameters, whereas McLachlan's principle requires 
that the error of many-body evolution must be minimized. On the other hand, 
the TDVP, as stated, requires stationarity of a given action. The three principles 
thus support quite different doctrines. 

Applying either the TDVP or the Dirac-Frenkel's principle, the authors in Q 
have proposed a method they called MCTDHB (Multi-Configurational Time- 
Dependent Hartrec method for Bosons). This approach has, for example, pro- 
vided tools for the description of the fragmentation of bosonic many-body states 
lol 11, 12, 13]. We will describe below in detail that, besides its many beneficial 



properties, the MCTDHB method is incomplete in certain situations. Specif- 
ically, when the single-particle density matrix (SPDM) becomes singular, i.e. 
noninvertible, the method fails. As a consequence, MCTDHB does not pro- 
vide a way to propagate, for example, a single condensate into a fragmented 
condensate many-body state. Although MCTDHB provides an important tool 
to describe the many-body physics of interacting bosons, the method therefore 
lacks the possibility to directly connect the phenomena of condensation and 
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fragmentation. 

Here, critically examining Dirac-Frenkel's principle and the TDVP, and 
adopting alternatively McLachlan's principle for truncated many-body evolu- 
tion, we improve on the previous multi-configurational Hartrcc methods, and 
solve the singularity problem of a noninvertible SPDM. In the process, we will 
also validate the resulting equations of MCTDHB in a different manner, how- 
ever additionally offering a straightforward handling of the exceptional evolution 
points related to the singularity of the SPDM. 

2. Variational Principles 

Let us now discuss the possible variational principles in more detail. We 
are aiming at finding an approximate solution of the many-body Schrodingcr 
equation when the state |$) is restricted. McLachlan's principle (jj, which was 
presented in 1963 as a new version of Frenkel's principle, requires the mini- 
mization of the error or remainder of this approximate solution from the exact 
evolution. The time evolution of any state is dictated by Schrodinger's equa- 
tion, idt\&) = H\Q). In other words, the evolution of state is determined by the 
Hamiltonian at any moment. But to make the state |4>) manipulable, we are 
generally forced to restrict or confine the state |$) into some simple and com- 
putationally feasible forms. With the state 1$) in restricted form, [idt — H]\$) 
cannot be exactly zero in general. Therefore, McLachlan's principle aims at 
finding the approximate solution which minimizes the positive semidefinitc er- 
ror measure ($| [idt— H] ''[idt — H ]|$). The details of the corresponding procedure 
will be rephrased in section 12.21 after introducing a concrete way to restrict the 
many-body state in a computationally feasible form. 

Hence it is guaranteed that the equation of motion obtained from McLach- 
lan's principle follows the exact evolution most similarily under given con- 
straints. The most appealing feature of McLachlan's principle is thus that it 
is a quite intuitive principle. Since it offers the possibility to evaluate the er- 
ror directly, we can intermediately, monitoring the error, increase the number 
of orbitals in the truncated held operator expansion, i.e. truncate the state 
less, to assure accuracy of the result. Alternatively, to save computational costs 
and time, the number of orbitals can also be decreased intermediately, i.e. by 
truncating the state more, in particular in cases where decreasing the num- 
ber of orbitals does not affect the accuracy of result significantly. Our scheme, 
described in detail below, in which McLachlan's principle is applied, therefore 
offers the opportunity to dynamically adjust the truncation of field operators 
or the state itself properly during computational time evolution, since we can 
monitor the error. As a particular benefit, a convergence test, mandatory for 
MCTDHB, is unneccessary, as we can directly obtain an error which indicates 
automatically how well our approach describes the interacting system of bosons. 

On the other hand, the TDVP which was carried out in MCTDHB H re- 
quires stationarity of action SS = 0. This does not necessarily mean an ex- 
trcmum (minimum or maximum) of the action. Though stationary points in- 
clude local extrema, the principle practically imposes only stationarity of the 
action: The equation of motion comes from stationary point of the action, which 
is not even necessarily a local minimum or maximum. 
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In many-body quantum mechanics, the action is written in terms of an 
expectation value of an operator-valued functional: 
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dt ($| 



id t ] + [id t ] f ] | 



(1) 



which is in quantum mechanical correspondence to the classical Lagrangian 
action. Here, Vt rap (x, t) is the (in general time-dependent) scalar trap potential 
confining the atoms, V{x a ,xp) is the two-body interaction potential, and m 
the mass of bosons. This is a real- valued functional of the many-body state 



|$) [l5(. Sometimes the action is simply expressed as S — f dt ($>\H — £<9t|$). 
Variationally changing the state |$) and the temporal change of the state t?t|$), 
we find the evolution of the state around the stationary action point. 

But here a problem occurs. When the form of the state is restricted or 
truncated for computational reasons in the sense that the state resides in a sub- 
Hilbert space, it is questionable whether the equation of motion obtained in the 
sub-Hilbert space leads to an evolution most similarily to the exact one obtained 
with the unlimited full Hilbcrt space. Even though the equation of motion 
obtained variationally with a non-truncated state can give the correct many- 
body Schrodinger equation, we cannot rely on the correctness of the equation 
in case the state is truncated. The path of stationary action with limitations 
imposed on the path can in principle deviate far from the one obtained without 
any constraints on the path. 

As a simple example, when we restrict the state to have only an overall 
phase change, i.e. |$(£)) = e~ in *|$), the action becomes S = J dt((H) - O). 
Depending on the value of fi, the action can be positive or negative. Actually 
there is no upper bound and lower bound on this action. Even worse, for some 
types of constrained states, there can be no stationary point of the action at 
all. So when applying the TDVP, at least the convergence upon increasing 
the number of orbitals, i.e., loosening the constraints, must be tested for every 
specific problem, to ensure the validity of the results. This is because, in contrast 
with McLachlan's principle, there is no direct error indicator in the TDVP which 
controls the accuracy of the approximation. 

The earliest variational principle for the approximate solution of many-body 
dynamics is Dirac-Frenkel's principle Q , which requires 

(5*\H-idt\*} = 0, (2) 

where (5$ denotes possible variations of the many-body state $ with respect to 
the variational parameters. The equation is quite similar to the TDVP when 
the action is given by S = J dt(Q\H — idt\&). The difference and (possible) 
equivalences between Dirac-Frenkel's, McLachlan's and other variational prin- 
ciples have been extensively discussed in the past 0, [l4[ . In Q , it is concluded 
that if the relevant manifold, i. e. the sub-Hilbert space, can be parametrized by 
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pairs of complementary parameters, the above mentioned principles are equiva- 



lent. In [14(, it is insisted that both Dirac-Frenkel's and McLachlan's variational 
principles are equivalent if both 5$ and <5$* are possible independent variations. 
But the "equivalence" merely implies the same resulting equation under some 
given conditions, not the equivalence of the principles themselves. In addition, 
it appears doubtful that <5<f> and 5&* are possible independent variations since 
5&* is simply the complex conjugate of (5$. Furthermore, as explained in detail 
later, principles which result in the problem of a noninvertible SPDM, which 
was alluded to already in the above, lack some information in comparison to 
McLachlan's principle, which resolves this problem. 

In summary, comparing the three variational principles, McLachlan's prin- 
ciple appears to be most suitable for finding a truncated many-body dynamics 
which approximates the real dynamics of interacting bosons. Adopting McLach- 
lan's principle, the variationally optimal truncation of the many-body dynamics 
of interactiong bosons can be adaptively controlled with a monitored error. 

2.1. Truncating a many-body state 

The limited or restricted forms of the state |<£>) for the truncated many-body 
dynamics can in principle take any form. In the simplest case, assuming that 
the occupation numbers of bosons concentrate in one orbital for the whole time 
of evolution, we can treat the many-body state with one single-particle orbital. 
More generally, the state will reside in a sub-Hilbert space of a specific form. 
An easily extendible and flexible form of the limitation on the size of the Hilbcrt 
space is the multiconfigurational time-dependent Hartree wavefunction ansatz, 
in which the many-body state of bosons is described as a linear combination 
of permanents |n), with a finite number M of orthonormalized time-dependent 
single-particle orbitals. Increasing the number M of orbitals, we can easily ex- 
tend the time-dependent sub-Hilbert space M. (t) . The basic steps to be followed 
in the procedure are as follows. 

Firstly, the many-body Hamiltonian is given by 



H = dx ¥ 



(3) 

Ax a Axp¥{x a )¥{xp)V{x a ,x^{x & )^{x a ). 



With a complete set of basis orbitals, the field operators of creation and anni- 
hilation of particles are expressed by the expansions 

oo oo 

¥{x) = ^a\cf>*{x) and t) = ^ a* , t). (4) 

i=l i=l 

Using these, the Hamiltonian can be written 

H = ' + 2 Y Vijkia\a\a k ai (5) 

i,j=l i,j,k,l—l 



where the single-particle matrix elements are 



dx 4>*{x) 



V 2 



2m 



<t>j{x), (6) 
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while the two-body interaction is represented by 



Vijkl = yy dfidf 2 ^*(^a)</>j(^)l / (x Q ,X (3 )(/)fc(f / 3)(/);(x a ). (7) 

We abbreviate sometimes, for the sake of convenience, 

4j = 4jkl = & la]a k a h oJt fc = a\a]a k , a\ jk = a\aja h (8) 



and so on. Then, the Hamiltonian can be written in short-hand form as H = 

h+±V = EZ=i e « a ij + 3 Si~,*,l=i 

To find the many-body ground state \G), we have to minimize the energy 
expectation value Eq = (G\H\G). Without any restrictions on the state |G), 
the actual exact ground state will be found. However realistically, we cannot 
describe a state exactly as this will in general require infinitely many orbitals. 
So we confine the state into a space of finite dimension, i.e. using only finite 
number M of orbitals, which is computationally feasible, 

\G M ) = J2 W ( 9 ) 

The above is regarded as the truncation of the many-body bosonic state. Here 
|n G Ai) indicates a normalized Fock state or a permanent 

\n)^ [ai) ^ = |vac) with f> = iV, (10) 

which is a particle state of which the individual members are composed of n\ 
particles in §\{x) orbital, n-i particles in (j>2 (x) orbital, • • • , and um particles in 
4>m{x) orbital. These M orbitals must be orthonormalized to each other. 

dx tf(x)^(x) =6ij. (11) 

These orbitals compose the sub-Hilbert manifold spanned by M orthonormal 
orbitals, which will be denoted by Ai in our context. The operators of creation 
and annihilation on these are related to the field operators in position space by 
the inversion of (@|, 

dx ${x)(/)*(x) and aj = J dx &(x)<f>i{x). (12) 

Here, a sub-Hilbert space spanned by c i4>i{x) is to be chosen so as to de- 

scribe the ground state optimally. And the coefficients C,-j which give minimum 
of energy are to be determined, too. Then, \G M ) can be considered as opti- 
mal truncation of the actual ground state \G). The details on how to proceed 
concretely will follow below in section [4j 

In the time-evolving case, we change not only the coefficients Cft along time, 
but also the sub-Hilbert space M. changes with time. With a time-varying 
truncation of the many-body state, we can express the state as 

i$(t)>= Yl c n(mt), as) 

n£A4(t) 



G 



with time- varying orbitals and their conjugate creation operators 

a\{t) = { &x¥{x)<f>i{x,t). (14) 



This approach was also incorporated in MCTDHB [§] . The time differentiation 
of the state |3>(t)) then becomes 



id t Mt))= [{idtCn(t))\n;t) + C a (t)f2 j d£ (id t <t>i(x,t))¥ \S)ai\n;t) 

1—1 

(15) 



n6M(t) 

Expanding dt4>i{x,t) with a complete basis yields 



id t (j>i(x,t) = y^J ki (t)(t>k(x,t). (16) 



k=l 

where the matrix tk% for 1 < k < M indicates an inner rotation inside the 
sub-Hilbcrt space, whereas t^i for k > M changes sub-Hilbert space itself. In- 
tegrating both sides after multiplying with <j)* k (x, t) , we obtain 

tu{t) = Jdx <f>l(x,t)(id t <f>i(x,t)). (17) 

Then Eq. (|15p can be expressed in the alternative form 

oo M 

id t \$(t))= [(id t Cn(t))\n;t)+CrMY,T, tk - & l^ t ^ ( 18 ) 

n£M(t) k=l i=l 

These preliminaries will be used in the following sections, by providing the tools 
for describing the truncated state |$) optimally. 

2.2. Adapting the number of orbitals 

The governing equation of MCTDHB which comes from either Dirac-Frcnkcl's 
principle or TDVP implies that the SPDM must be always invertible, not only 
initially but also at any instant afterwards. We quote here the equation (26) in 
[911, which is 



Z=M+1 



d t (j)k{x 1 t)= ^ ei k + ^ (P) kt (ali pq ) v inp 



M 

'<t>l&t) (W) 



i,7l,p,q—l 



in our notation, where {p)^i represents the inverse of the SPDM (pij) = 
In it seems to be taken for granted that the SPDM is always invertible. 
When the SPDM becomes noninvertible, however, the MCTDHB method in 
fact suddenly fails. For example, when the whole bosons of our interest resides 
initially in single orbital, MCTDHB cannot propagate this pure condensate into 
any fragmented state. As a consequence, the method does not provide a way to 
find the form of the second orbital. 

To decrease the sub-Hilbert space or the number of orbitals is not an issue. 
We can simply eliminate those orbitals with an occupation number not of order 
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N. However, to increase the sub-Hilbert space dimension, that is the number of 
orbitals M, becomes difficult since an increase of the sub-Hilbert space and an 
additional orbital must be set up optimally. Noninvertibility of the SPDM can 
obviously also happen dynamically. Unoccupied orbitals or scarcely occupied 
orbitals thus cause a problem with the dynamics of MCTDHB. These prob- 
lems cannot be resolved by Dirac-Frenkel's principle or the TDVP. As McLach- 
lan's principle is based on requiring that the least error is acquired during time 
evolution, we may resolve this problem by finding an additional orbital which 
minimizes the error, as will be expounded in detail below. 

McLachlan's principle does not use the action, but the concept of error or 
remainder. We know the exact form of the many-body Schrodingcr equation. It 
is idt\$>) = H\G>) with the full Hamiltonian Eq. ([3]) in second quantization form. 
Since the exact calculation is too cumbersome, we can represent the state by 
the relatively simple multiconfigurational time-dependent Hartree form. The 
remainder from exact evolution in any case becomes 

[idt-H]\*). (20) 

In this expression, the left part idt\&) is an evolution of the state |$) in its 
limited, truncated form and the right part H represents the exact evolution. 
Here, the initial state, for either part, is specified in its truncated form, while 
the Hamiltonian H is not truncated. 

A quantitative measure of the instantaneous error is then 

($\[id t - [id t - H]\$) >0, (21) 

which is (by definition) positive semidefinite. We have to minimize this error 
by varying the truncated state in Eq. (|13|). 

3. Method of steepest constrained descent 

In this section, we describe the basic method which will be used to find 
both a minimum of energy and a minimum of instantaneous error in the fol- 
lowing sections. Since we use complex variables and complex functions, the 
well-known method of Lagrange multipliers will here be carefully extended to 
one with complex numbers. In addition, we introduce the method of steepest 
constrained descent, which propagates an initial trial estimate (the many-body 
state) by steepest descent under the given constraints toward an approximate 
minimum point. The method provides a simplified way to find the latter, and 
is introduced because a fully self-consistent solution is frequently forbiddingly 
difficult to find from the equations following from using the bare method of 
Lagrange multipliers. 

3.1. Method of Lagrange multipliers with complex numbers 

In mathematical optimization, i.e. when we want to find a minimum or max- 
imum of a continuous differentiable function f(x) subject to some constraints 
gi(x) = const, where the index I enumerates the constraints, the method of La- 
grange multipliers is widely used. Providing conditions for stationary points 
subject under the constraints, it provides a strategy for finding minima and 
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maxima. Comparing the obtained values at the stationary points obtained, we 
can determine which point is a minimum or maximum. 

With complex variables or complex functions, some caution is however nec- 
essary. As one complex variable Zk corresponds to two real variables afc + ibk, 
both real and imaginary parts of z k should satisfy the Lagrange formula, 



oa k , oak 



Ob k V ob k 



(22) 



where / and gi are given real-valued functions, 
express the two real equations above. 



df_ 

dzl 



i K 



o.Fm,. 2idb k ^ 1 



2 dak 
1/0/ 



l 77 



2\da k 



2da fe 



There is a compact way to 



2i db k 

-E A 



96, 



(23) 



Since the functions / and gi's are real- valued, the derivatives with respect to real 
variable a^'s and bkS are also real- valued functions. Hence the above complex 
equation requires that both real and imaginary parts would be zero indepen- 
dently This is not because Zk and z^, are possible independent variables. Zk 
and its conjugate z^ are definitely dependent variables. We use this just as a 
compact way of expression for two real equations, using partial derivatives with 
respect to the complex variable. Here, we need to be very cautious using the 
partial derivatives. When a function / is given with several variables, the par- 
tial derivative df jdz means a differentiation with respect to only that variable 
z even though others such as z*,a,b are dependent on z. We should also be 
careful with which variables the function is expressed. Commonly df/dz means 
df(z,z*)/dz and df/da means df(a,b)/da. Only when there is no possible 
source of confusion, we drop these explicit indications for the sake of conve- 
nience. 

The formulas change if the constraints gi are given in complex form. A single 
complex constraint corresponds to two real constraints, 



9l 



(R) , ■ (I) , 

g\ + ig\ = const. 



(24) 



Then, with two groups of real Lagrange multipliers for the real part and imagi- 
nary parts of the complex constraints, 



df. = 

dx k ^ 

= E 



dR) dg\ H) | x {i)dg\ I] 



dx k 





dgi _ 




A l) 


'dgi 






2 


-dx k 


dx k - 


+ 2i 


-dx k 


dx k - 







dgi , 

dx k 




A«- 


dgn 


2i - 


- 2 


2i - 


dx k 



(25) 
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(26) 



Complex constraints come with complex Lagrange multipliers. Then, the com- 
plex conjugate form of the constraints must be added. When both variables and 
constraints arc in complex form, the equations will be represented by 



Case by case, we can choose the simpler form of the two equivalent complex 
equations (|2T|) and (|2"5)l above. 

3.2. Approaching along steepest constrained descent 

Even when we have found the equations for extreme or stationary points, 
we usually cannot find those points directly from the equations except for some 
simple cases. What we conventionally can do at most is, in general, to give some 
approximate estimates of an extremal point. These estimates would not satisfy 
the stationarity condition exactly. Therefore, to find the actual accurate point 
of extremum which fulfills the stationarity condition found by the method of 
Lagrange multipliers, we have to propagate the initial trial point of departure, 
which is in a truncated many-body state space in our context, to the actual 
minimum or maximum point. 

Allocating a pseudo time r along the path on which the point in the config- 
uration space (the many-body state in our context) moves, we can consider the 
rate of change 



with respect to the pseudo time r. The rate of change -f- indicates how rapidly 
the value of / changes at x(t) along the path toward t 2 . The Lagrange formula 
Eq. (|22|) comes from the condition that 4*- = for any variation of the Xk, i.e. 
any ^p-, under given constraints. In non-stationary points, the path which gives 
the smallest (biggest) value of $f- will be called the steepest descent (ascent). 
However, without any limit on the velocity of the variable changes t 2 , the rate 
of change 4^ would be also unlimited. In addition, any given points along the 
path must still satisfy the constraints gi(x(r)) = const, as we are trying to find 
the minimum or maximum of f(x) under the given constraints. So minimizing 
(maximizing) the rate of change df /dr Eq. (|29p within the constrained path, i.e. 




(27) 



or equivalcntly 




(28) 




(29) 




(30) 
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and under a closed topology which can be, for example, taken to be a hyper- 
sphere 

const, (31) 



(dx k V- 



we can propagate the point along the steepest constrained descent (ascent). 
Here, what we are varying is the speed of variable change, i.e. the set {^r}- The 
closed topology which we have chosen limits the overall rate of change (speed), 
of the variables |4p|, so that the rate of change ^ does not diverge. Without 
a closed topology which limits the change of variables, we would not be able 
in general to find the steepest path since it possibly diverges with the variable 
changes under an open topology. To be specific, in our context, the function 
/ to be minimized corresponds to the energy expectation value E g m (or the 
instantaneous error (\id t — Hy\idt — H])), the constraints gi to J^n^n^n = 1 
and J dx <f>*(x)(f>j(x) = Sij for i,j < M (or time derivatives of them); therefore 
I goes from 1 to (1 + A/ 2 ), and the variables Xk correspond to the coefficients 
Cfi and the orbitals <f>k (x) (respectively their time derivatives) . 

As we minimize (maximize) the rate of change df/dr Eg. (f2U]) under the 
closed topology with respect to the rate of change, ^p-, we use again the method 
of Lagrange multipliers which gives the condition of the stationary 4*- in the 
form 

!^-E^-^ = 0, (32) 
dx k ^ dx k dr 

for all possible k. The simple linearity of Eq. (f2THl3TI)) on provides directly 
the result 



dx k 



df _y- A 

dx k 1 dx k 



(33) 



from Eq. (|32j) . Here we transformed the multiplier A into a new convenient 
expression ^A(r) = 2 \(t) ■ ^he undetermined Lagrange multipliers A/(r) and 
A(r) are then determined by the constraints Eq. (j30|) and Eq. (|3Tj) . Inserting 
Eq. (122|) into Eq. O, 



-r- A/ \r d 9m 



k 



df y^ A dgi 
dx k 4^ ' dx k 







E g/ dg m _ x \p 9gi dg, 
fin-, Fin-, ~ 2^, 1 



(34) 



dx k dx k ^ , dx k dx k 

k 



Introducing the inverse of the symmetric matrix A\ m = gjr-T^FS A/ is 
determined to be 

df dg m ! 

= f-' ~dx~k ~dx~k ml (35) 

Though A(t) is also to be determined to satisfy Eq. (f3"Tj). the 'const' on the 
right side of Eq. (f3"Tj) is just a constant whose value is unimportant. Unless we 
fix the size of the topological constraint 'const', A(r) can be whatever function 
of t. But since the rate of change ^ should be negative (positive) to approach 
a minimum (maximum) of f(x), we choose A(t) to be always positive definite, 
i.e. A(t) > for all r. Furthermore, we take a minus (plus) in front of the 
right-hand side of Eq. (|33[) to approach a minimum (maximum) of f(x). 
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The point in configuration space will eventually arrive on a neighboring local 
minimum or maximum along the path 



f T dx 

'(t = T) = x(t = 0)+ / dr ^(r), 
Jo dT 



(36) 



where Xk{r = 0) is the initially estimated point and ^f(r) represents the direc- 
tion of the steepest constrained descent (ascent) under a given topology. The 
path will strongly depend on the choice of topology. To save calculational time, 
we can fix some variables alternatively. In other words, setting dxi/dr = for 
some set {/} within some pseudo-time r span, we change only the other vari- 
ables along the steepest constrained descent. Then the hypersphere topological 
constraint Eq. (|3Tj) on the velocity of the variable, changes into 

dxi\ x -» /dxk\ 2 

k=£{l} 

It approaches in any case a point which has less (more) numerical value of 
f({xk}) than the current point, following a stepwise path. As a simple example, 
when /(a, b) is a function of two variables, we can change a and b separately to 
approach a neighboring minimum. When §^ > 0, we decrease a until = 0. 
After then, testing ^ at that point, we propagate b until §{ = 0. Successively 
and iteratively we change a and b until both f~ = and ^ = are met. 
Along this stepwise path, the point propagates to the neighboring minimum. 
By this means, which is a separation of the change of variables, we can obtain 
an efficient and thus fast procedure for a concrete computation. 

As the steepest descent (ascent) does not gurantee the shortest path to the 
minimum (maximum), a flexible and adaptive choice of topology is crucial to 
save time and effort. However, the hypersphere topology and separation of the 
change of variables will be sufficient in most practical cases. 

In summary, the method of steepest constrained descent (ascent) guarantees 
that the configuration space point approaches a neighboring local minimum 
(maximum) of f(x), by propagating the point toward 



dx k , 
~*F = TA(T) 



df ^ x dgi 

dx k 1 dx k 



(38) 



for hypersphere topology. Here A/(r) is given by Eq. (|35p at each point, and A(r) 
is any arbitrary positive function of r, which would be chosen in any convenient 
way during the actual computation. To approach a minimum (maximum) of 
f(x), we take minus (plus) in front of the right hand side of Eq. (l38l) . 

With complex variables and complex constraints, using Eq. (|27l) . this will be 
simplified as 



dz k da k .db k , 

— = — h i—— = tA(t 

dr dr dr 



(39) 



We have ignored that A(r) should be doubled on the right-hand side since A(r) 
can be arbitrarily chosen, i.e. we put 2A(r) — > A(r). 
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4. Minimizing the energy of the truncated many-body state 



4-.1. Variational method with Lagrange multipliers 

Before investigating the dynamics, let us demonstrate the time-independent 
scheme first. To find the many-body ground state, we have to minimize the 
energy expectation value 



E GM = (G M \H\G M ) = [ J2 HQ]£[ E C «l 



HeM 



(40) 



by varying the coefficients C^'s and the set of M orbitals subject to the 

(1 + M 2 ) constraints 



y] C~C.it = 1 (1 constraint) 



(41) 



and 



dx 4>i(x)<j>j(x) = 8ij iori,j<M (M constraints). (42) 

The variation with respect to the expansion coefficients gives 
8E g m 



dC% 



= (n\H E C,n\m) = XCn = E g mCh, 



(43) 



where ft <G M. We have a real functional E g m Eq. ([4T)| to be minimized, and one 
real equation of constraint Eq. (|41j) . with ^^l-iy. com pl cx variables Cn when 
we fix the total number of bosons N . The undetermined Lagrange multiplier 
A which has to be real is determined to be E g m with the help of constraint 
Eq.gU). 

Using the properties e£ = e^, V ljk i = Vjiik and V* kl = V ikjl = V kHj , we can 
express the above equation explicitly as 



niCn 



H E C^\rh) — E g mCh 

E e » + ^ V uu(ni - 1) + ^ E ( Vlklk + V lkkl)t 

I L k 

E ' £ y + V m( n l - !) + V ljjj n j + ^'{Vi k kj + Vi k]k )n k (n 3 - + \)n x C^ 

j,l L k J 

^ 2 H' V ^o V /K+2)(n J + l)(n i -l)n ; C^ J 

i- 2 S ' v // ^ + 1 )(^ + 1 )( ni - 1 ) niC '^ 
h 2 X ' ^ y / («j + 2 )K + 1 )«fc^c^ 



(44) 
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where the primed summation is performed such that different indices can 
have only different values. For example, summations over two and three indices 
are Ej,z ' = Ej Ej^-j ' = E* E j¥ i Ei^i or j. and E fe ' insid <3 a bracket 

means Efc^(the other indices)- Though Eq. gU) is just an eigenvalue equation with 
fixed matrix components, the terms and Vijki are matrix elements depending 
on the orbitals, which are to be determined as follows. 

For variation with respect to the orbitals {<f>k}, functional differentiation is 
used. As the permanent |n) is constructed from repeatedly applying creation 
operators of particles in the M orbitals, it can be regarded to be given by 
multiple integrations over the orbitals <pj.(x), 



(45) 



As each functional differentiation contributes to the result, this will be counted 



dx a 4>i{x a )^ (x a ) / dx^j(xp)^(xp) ■ ■ ■ | vac) 



by a factor at , which results in 



d(<S>\ 



$|a[,\E , (x). But as the full hamiltonian 

H, which is given in field form by Eq.(|3]), is independent on the M orbitals 
chosen, the functional differentiation of H with respect to the orbitals {<fik} 
gives zero, g®! 1 ^ = 0. Then functional variation of E g m Eq. (|4T)|) with respect 
to the orbitals {<fik}, combined with the functional constraints Eq. (|4"!?|) . leads to 



dE, 



d(j>* k {x) 



.E< 



M 



3=1 



(46) 



where Ajk = Ajy is Hermitian matrix. Here the method of Lagrange multipliers 
with complex functional variables is used. Integrating each side over space after 
multiplication with (j>*(x) yields 



aiaiH 



[ E C *\ H < 



for / < M 
for / > M. 



(47) 



Here, using Eq. (|43p and the property that m l h € M. when k,l < M and m € Ai, 
the undetermined set of Lagrange multipliers X^i becomes related to Eqm by 
Aft/ = E g m{o\,i) for k,l < M. Using Eq. and the bosonic commutation 
relations between field operators, i.e. [&{x), &(x')] = 5{x — a?), [^(x),a\,] = 
<f>k(x), and so on, Eq. (|4"6")) {a^f(x)H) = Ej=i ^kj4>j(x) is explicitly expressed 

(&{M(x)) 



as 



EO 



V 2 

- 2^ + V ^vW 



4> 3 (x) 



E < 



df ^(x , )y(f,f)^ g (f)0 J (f) 
Since we can eliminate the annihilations above M, we obtain 

E (4L)^(^) = Eh. 



M 

^ \ k j<t>j(x). 



(48) 



M 

E< 



a/ 



i=i 



M M 

(4L>^(*) = E 

p,g,3=l 3=1 
A/ 



- (at Fa 



(49) 



G m {a\.hj 



(alHaj) 



(f>j(x) 



3=1 
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where the new undetermined Lagrange multipliers Xkj = Xkj — (a^Haj) which 
satisfy A£ • = Xjk are introduced. We abbreviated, for the sake of convenience, 
two single-particle and interaction operators defined by 



h(j)j (x) 

and 



2 



V 



h(x) (50) 



V pq <t>j(x) = J Ax' <^(f )V{x,g (51) 

^.U. Approaching the many-body state 

To find the tentative ground state |G M ), we have to find the coefficient C^'s 
and the complex orbital functions (f>k(x) satisfying Eg. (|4"4")l and P^|) . Addi- 
tionally, the solutions must satisfy all (1 + M 2 ) constraints Eq. ([4"Tll4"2"j) . The 
number of real values which we should find is 2 yj7^rm for the set of the Cs, 
2M real functions for the </>fc(x), and (1 + M 2 ) for the undetermined Lagrange 
multiplier E g m and the Xkj. The number of given real equations therefore is 
2 ^mW^TJT for Ec l- B. 2M real functional equations for Eq. (g5|, and (1 + M 2 ) 
for Eq. (|41U42j) . Apart from the large number of variables, what makes matters 
even more complicated is that the equations (|4Tll42ll44ll49)) are coupled together. 
To find a self-consistent solution is therefore obviously a very difficult problem. 



In ref. [17|, the authors started from an initial guess, then iteratively, with 
a convergence check, they obtained a solution. A few years later in ref. @, 
applying the Wick rotation it — > r on the equations of motion, they introduced 
the so-called imaginary time propagation. They claimed that this reduces any 
arbitrary initial many-body state after a sufficient time of propagation to the 
ground state. The imaginary time evolution i§^\^) = H\&) =>■ -^|$) = H\$) 

implies that e~ lHt => e~ Hr \&). As r goes to infinity, this seems to indicate that 
only the ground state survives and the excited states would no longer contribute. 
The contribution of the excited states decays exponentially according to a factor 
which is proportional to energy difference from the ground state and r. 

Here, we present alternatively the method of steepest constrained descent 
described in section [3] This method of a approaching a minimum (maximum) 
from an initial trial state will guarantee that any given state is propagated to 
the neighboring lowest (highest) value point along the steepest path for given 
variational parameters. Though it propagates a state only to a local neighboring 
minimum (maximum) point, we can find in many cases the global minimum, 
or the ground state, from a well-chosen initial state, and with in addition well- 
chosen variational parameters. Although this does not deliver the state along 
the shortest path, the very large number of degrees of freedom on the choice 
of variational parameters or the sequential processing (separation) of variations 
can compensate it in many cases. The flexibility of this method will therefore 
be beneficial for finding the ground state. 

Let us see the process in more detail. The first step is to find (by educated 
guess) an appropriate initial state, specifying the coefficients and the M 
orbitals which we believe are appropriate to approximately describe the ground 
state of a given system. From this initial guess, the state is propagated as 
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follows. For the expansion coefficient Cfj's, using Eq. (|39l) . 



(It 



-A c (t) 



|A^Cfi|7n) - AC*, 



(52) 



where Aq(t) is any arbitrary positive function of r which is introduced to satisfy 
^2fi ^T^§T = const at a certain given instant r and therefore can be chosen 
in a convenient way to save time and calculation costs. Since it must satisfy 
En CtC n = 1, i.e. £ s + d -§fC n ] = 0, where A = (H). 

As another option, we can use a polar representation for Cfj. Representing 
the complex variable Cs by an Eulcr representation with a radius £s and an 
angle On gives = £se ie ". Then the constraint Eq. (|4Tj) becomes 2n£l = 1) 
restricting only the radial component of the complex variables Cf{. Since we 
do not have to confine the variable change into the specific form [(^-) + 
£| (^f) 2 ] = const, separating the two variable sets can be much more efficient in 



this case. That is, we set = and V- 

' dr t-~in 



const for some given r time 
spans which would be chosen in computationally favorable way, and -jA = 

and (^dr~) 2 = const for the other t spans. Then the method of steepest 
constrained descent gives 

-7- = 0, -7- = -Afl(r)9(foe 



(53) 



for some given r spans, and 



d9n 
dr 



= 0, 



dr 



= -Ac(r) 



1^)) - A& 



(54) 



for the other following T-intervals. Here 5ft means real part of complex num- 
ber and 9 means imaginary part of complex number. Using = 1j an d 
therefore Yln^^d = 0, then A becomes A = (H) too. Dealing with ^ and 
#s separately, we propagate the two variable sets successively and iteratively 
until convergence is achieved. With any sequence, they will finally approach 
the minimum. 

For the orbitals 4>k (x) , the method of steepest constrained descent gives 



d(j>k{x) 



dr 



M 



= -A 0fc (r) (a{*(x)H) - Yl 



-A 0fc (r) 



M 



j'=i 

M 

E ( 



(55) 



>V m (j>j(x) - ^2 A fc: ,^(x) 



If we propagate the orbitals separately one after another, i.e. only the fc-th or- 
bital changes within a certain period, the constraint becomes J dx (f>* (x) 



4>k(S) 

dr 



for I ^ k and / dx {4>* k {x)^f- + ^f-(j) k (x)) = 0. Then the undetermined 
Lagrange multipliers becomes 



M 



M 



(56) 



P,9J = 1 
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and 



M 



M 



l kpqj 



l"IJ 



(57) 



Separating the propagation of the orbitals, i.e. propagating the orbitals one after 
another independently, makes the process much easier and more comprehensible. 

For the numerical implementation by computers, a discrete variable repre- 
sentation technique (commonly abbreviated DVR in the literature) will be used 
to represent the continuous orbital functions with a finite number of discrete 
variables fl6j . Based on well-known orthonormal functions such as, e.g., the har- 
monic oscillator eigenfunctions, sinusoidal, Bessel, and Legendre polynomials, 
we can express any continuous function as 



1=1 



(58) 



This is mathematically founded on the completeness of eigenfunctions and Sturm- 
Liouville theory. With an appropriate choice of the eigenfunction set, several 
hundreds of eigenfunctions will be sufficient to describe the continuous functions 

4> k {x) [2 EES [II . 

Then the change of orbitals will be expressed as 



^7 = - A ^ (r) 



M 

E<4 

0=1 



kj/ e lj 



M 

E 



M 



stt 

l kpqj/ v lpqj 



3 = 1 



(59) 



With these successive and iterative propagation steps, we can approach the 
ultimate form of the ground state \G M ) in the sub-Hilbcrt space. 



5. Control of truncated many-body evolution 



5.1. Evaluating the error of truncated many-body evolution 

An instantaneous error is expressed by Eq. (|21j) . Minimizing this error with 
a state change under the truncation Eq. (|13|) gives us the appropriate truncated 
many-body evolution. This offers, as a major benefit of the present approach, 
a quantitative value, the error, which indicates how accurately the truncated 
evolution describes the exact one. 

Explicitly expressing the error, we have 



®\[H-id t ]\H-id t ]\3>) 

n\C*,H + (fiUidtCi) 



= E 
x E 

meM(t) 



M 

E' 



n\Cta] 



M . 

HC A \rh) - (id t Cm)\m) ~ E / dx" {id t <f> 3 (x' ,t))¥ (x')a 3 C^\ 



(60) 
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We minimize this instantaneous error, varying the complex variables (ptCfi) 
and {pt4>i{x ',£)) subject to the (1 + M 2 ) constraints 



£Q(t)CW(i) = l 



E 



{d t Cm)Cn{t) + Cm{dtCn(t)) 



(61) 



and from the orthonormality condition 



dx<j>i(x,t)(l>j(x,t) = 5ij 



d.r 



# (a?, tjjfa (x, t) + 0* (x, t) (dtfa (x, t)) 



(62) 



The Eq. (|5T|) requires probability conservation of the state itself and Eq. (|6"2"j) re- 
quires conservation of orthonormality of orbitals. Using the expression Eq. (|17[) . 
Eq. (|62[l can be expressed as the hermiticity condition tji = tij . 
Variation with respect to dt CS leads to 



{n\i[H - id t }\§) = X(t)C n 

resulting in 

oo M 
i=l j=l 

As the constraint Eq. (|61|) enforces A(i) = 0, 



(63) 



(64) 



(65) 



(n| [H-idt]\$) =0, 
and therefore the time evolution of the expansion coefficients takes the form 

id t Cft = (n\[H-t]\<f>) (66) 

where n £ M(t) and i = J^i j UjO>\j. 

Variation with respect to 9 t t) gives 



($|za£*(2) [H - id t ] |$) = A «<M^ t). 



(67) 



After multiplying both sides with </> ; * (x, t) and integrating, one finds 



{ma\ai\H - id t ]m = 



X H (t) itl<M 
if I > M. 



(68) 



Since {n\a) kl belongs to M(t) whenever n £ M{t) and k, I < M, all Xki(t) here 
become zero too with the help of Eq. (|65"T). So for any k and Z, 



(*\&l l [H-id t ]\$)=0. 



(69) 
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Expanding the above equation for I > M 

M 



(70) 



Since the SPDM can be noninvertible, we reduce the density matrix by elimi- 
nating unoccupied orbitals in which the eigenvalues of the SPDM becomes zero 
(relative to O(N) in the limit N — > oo) after diagonalizing the SPDM. In the 
process of diagonalizing the SPDM, the orbitals are unitarily transformed so 
that the essentially unoccupied orbitals can be found and eliminated. Introduc- 
ing the inverse of this reduced density matrix (pki) = = wnen the 
O(N) occupied orbitals exist up to the Mith orbital (Yli=i(p)ki~ (Pij) = $kj) an d 
using the completeness relation ' ■> i (a?, t) = 8{x' — x), the evolution 
equation of the orbitals acquires the form 



A I 



oo Mi 



1=1 l=M+l i=l 

AI oo All r Mi All 

Y,tikMx,t)+ E Z>>«< E e <« a t«+ E V 

1=1 l—M+1 t=l L„ =1 n,p,q=l 

M oo r Afi 

t" ) 



E^fc0z(^,i)+ E 



;=A/+l 



i,n,p,q—l 



a 11 



<t>i(x,t) 



)0i(x, t) 



(71) 



for k < M x . 

Here, we divide the evolution of the orbitals into two parts. We call the 
left / < M part of Eq. ((TT|) inner rotation and the right I > M part of Eq. (jTTj) 
rotation toward the outside the sub-Hilbert space. Since the sub-Hilbert space 
spanned by M orbitals, J2i=i t), does not change under inner rotation, 

we realize that only a rotation toward the outside deforms the sub-Hilbert space. 
For the evolution inside the sub-Hilbert space, i.e. inner rotation, we simply use 
tij defined in ([TCI) , which can be any Hermitian matrix. Using Eq. (|65H69llTTj) , 
the expression for the error Eq. (|60p can be strongly simplified: 



oo 

[A -w t ]\*) = {<i>\[- J2 ^«ag w ][#-ift]|*} 

1 M± oo 

= o E E v ijkl ($Hh[H-id t }m 



i,j=l k,j=l 

All oo 

E E V wVlspq(aljk n )(p)n 

i : j,k,n : r,s : p,q=l l=M-\-l 

All AI oo ^ Mi 

' E E E V ijki v ls P g(al f jk at pg ) + - E E V m v ik. 

i,j,P,q=l k,s=l l=M+l i,j,p,q=l k,l=M+l 



S nr (®rspq) 



'q\®ijpq) 



(72) 



The above equation represents our main result, rendering the error of many- 
body quantum evolution upon truncation systematically computable. As clearly 
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seen, the error stems entirely from interactions. In other words, the error does 
not depend on the choice of Uj and even on the single particle energy matrix 
for any given truncated initial state, provided the evolution of the state, dt |$) , 
is optimally taken as in Eq. (p5l [551171]) . 

5.2. Determining the number of orbitals dynamically 

When the error becomes large or, alternatively, when we aim at describ- 
ing the system more precisely, we have to increase the number of orbitals Mi 
into M. The (M — Mi) additional orbitals then can be determined by varia- 
tion of the error with respect to <j) u {x,t) where Mi < u < M, and subject to 
the ortho normalization condition J dx </>Jj(x, t)<p v (x, t) = S uv . The method of 
Lagrange multipliers for functional variables gives the stationary condition 

Mi 

Y (V sp (f) q (x, t)) (a|t fer[ ) (p)-} ( 

i,i,k,n,r,s,p.q—l 

73 

Mi M M y ' 

- Y ]C V ijku{V sp 4> q {x,t)){a\l k a\ m ) = ^2,Huv<i>v{x,t). 

i,j,p,q=l k,s—l v—1 

While this is an equation for the additional orbitals, the additional orbitals are 
difficult to obtain directly from the above equation. As mentioned already in 
section IB~2l the best we can do is to guess some orbitals that will approximately 
satisfy the above equation. 

Therefore, we use the method of steepest constrained descent again. From 
an initial trial orbital, we propagate the orbital toward 

d(/)u(x) , / 



(It 



Y V Hku (V sp <f) q (x, t)) (aljkn)(P)nr 

i.j,k,n,r,s,p.q—l 

(74) 

Mi M M -| 

- Y Y Vijku (V sp <j) q (x, t)) {aljk&tpq) - Y ^nv<Pv{x, t) 
i>JiPi<l=l k,s=l v—1 

so that the error become minimized with this additional orbital. If we prop- 
agate orbitals separately and iteratively one after another, i.e. only the itth 
orbital changes along r, the constraint becomes / dx </>* (x) = for v ^ u 

and J dx (0* (x) + ^^-(j)u(x)) = 0. Then the undetermined Lagrange 
multipliers become 



Mi 



(J-uv — Yj ^ijkuV vsp q(aljf. n ) (p) 



Va" ) 

nr x^rspql 



i,j,k,n,r,s,p,q=l . . 

Mi M ^ ' 

— Y2 VijkuVvspq ( a ijk®'spq)' 

i,j,p,q—l k,s—l 

where the orbital index rangea are constrained to be Mi < u < M and 1 < v < 
M. 

Since the inner rotation can be arbitrarily chosen, we can take the unit 
matrix, iy = e^, which renders the result in a simple form. The evolution 
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Eq. (|66p for the expansion coefficients becomes 

M Mi 

id t C H = (n\ - J2 E V w4w l*> • (76) 

i,j=l k,l=l 

This implies the desired property that, when the interaction is turned off, the 
coefficients C r t do not change at all. The Schrodinger equation for the orbitals, 
Eq. (|7T1l . becomes 

oo Mi 

id t 4> k {x,t) = h<t> k {x,t)+ J2 E V inpa(p)u(att Pq )<t>l(Z,t) (77) 

;=M+1 i,n,p,q=l 

for fc < Mi. The projection toward the outside of the sub-Hilbert space M. 
takes place only on the interaction term. Thus the option of taking t\j = e,j 
for the inner rotation shows the effect of the interaction term in this explicit 
manner. 

For Mi < k < M, the evolution of the additional orbitals can be chosen 
in any convenient way, since the time derivatives of the additional orbitals do 
not occur in the error measure. The only constraint is the inner rotation. As 
we have chosen tik( = J dx(f>^(x, t)idt<f>k{x, t)) for Mi < I < M and k < Mi 
to be eik, tki should be t* lk = e\ k = e k i to satisfy the orthonormality constraint 
Eq. (j62l) . Taking the rotation toward the outside of the sub-Hilbert space to 
be also equal to the single-particle energy matrix, e k i for k > M, the evolution 
of the additional orbital <fii (x, t) for Mi < I < M is determined by the simple 
equation 

id t (pi(x,t) = hcf>i(x,t), (78) 

while the additional orbitals are found from Eq. (|74l) . 

We, finally, emphasize again that the key difference from the MCTDHB 
approach is that additional macroscopically occupied orbitals during dynamical 
evolution can be found in our approach. By this means, we can handle the 
exceptional case when the SPDM is not invertible, and increase the number of 
orbitals under any given circumstances and boundary conditions. 

6. Summary 

Using McLachlan's principle, which is physically intuitive and simple, and 
the methods of Lagrange multiplication and steepest constrained descent, we 
have found a computationally feasible way to describe the many-body evolution 
of bosons in a rigorously controlled manner. Writing the many-body state in 
Hartree form and limiting the size of the Hilbert space by truncating into a 
finite number of macroscopically occupied field operator modes, the error from 
the exact evolution can be minimized self-consistently. This gives a variationally 
optimized solution to the evolution of the truncated many-body state. 

We have demonstrated that without two-body interactions, our scheme pos- 
sesses the desired property that the evolution of the many-body state can be 
exactly described with zero error, cf. Eq. ([72)1 . When interactions are turned 
on, the error accumulates during time evolution. Employing our method, we 
can evolve the truncated many-body state in an optimized way. Monitoring the 
error simultaneously, we can increase the accuracy of the evolution by increasing 
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the number of orbitals in a self-consistent way. By adaptively changing the num- 
ber of orbitals based on the instantaneous error measure, we can automatically 
ensure the validity of the result for the many-body state. 

We conclude by a brief summary of our approach applied to the well-known 
and ubiquitous Gross-Pitaevskii equation. We start by evolving the initial trial 
state |$) along the M = 1 version of Eq. (|76[I77|) and simultaneously check 
whether the error Eq. (|72p remains small or not. Monitoring the error Eq. (|72[). 
we can determine under which conditions the Gross-Pitaevskii equation loses its 
validity. When this happens, the error becoming large, we increase the number 
of orbitals to M = 2 (thus, here, Mi = 1 and M = 2). The additional second 
orbital is found by the method of steepest constrained descent, using Eq. ([74)). 
Then the subsequent evolution of the quantum many-body state follows the 
M = 2 version of Eq. ([751177)) , while the evolution of the (initially singular) 
second orbital follows Eq. ([75)1 . We monitor the error Eq. ([72")) again, checking 
that the error is decreased to a sufficient degree. In a self-consistent manner 
the processes described are carried out successively until we obtain a prescribed 
accuracy. 
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